home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / tests / matode.tst < prev    next >
Text File  |  1999-09-16  |  9KB  |  308 lines
  1. Leps=1.e-6;
  2. //     dy1/dt = -.04*y1 + 1.e4*y2*y3
  3. //     dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
  4. //     dy3/dt = 3.e7*y2**2
  5. // on the interval from t = 0.0 to t = 4.e10, with initial conditions
  6. // y1 = 1.0, y2 = y3 = 0.  the problem is stiff.
  7. //     definition of rhs
  8. deff('[ydot]=f(t,y)',...
  9.    'f1=-0.04*y(1)+1e4*y(2)*y(3),...
  10.     f3=3e7*y(2)*y(2),...
  11.     ydot=[f1;-f1-f3;f3]')
  12. //     jacobian
  13. deff('[jac]=j(t,y)',...
  14. 'jac(1,1)=-0.04;jac(1,2)=1.e4*y(3);jac(1,3)=1.e4*y(2),...
  15.  jac(3,1)=0;jac(3,2)=6.e7*y(2);jac(3,3)=0;...
  16.  jac(2,1)=-jac(1,1);jac(2,2)=-jac(1,2)-jac(3,2);jac(2,3)=-jac(1,3);')
  17. //    
  18. y0=[1;0;0];t0=0;t1=[0.4,4];nt=size(t1,'*');
  19. //    solution 
  20. yref=[0.9851721 0.9055180;0.0000339 0.0000224;0.0147940 0.0944596];
  21. //
  22. //  1. fortran called by fydot, without jacobian
  23. y1=ode(y0,t0,t1,'fex');
  24. if maxi(y1-yref) > Leps then pause,end
  25. //  2. fortran called by fydot, type given (stiff), no jacobian
  26. y2=ode('stiff',y0,t0,t1,'fex');
  27. if maxi(y2-yref) > Leps then pause,end
  28. //  3. fortran called by fydot , fjac, type given
  29. y3=ode('stiff',y0,t0,t1,'fex','jex');
  30. if maxi(y3-yref) > Leps then pause,end
  31. //   hot restart
  32. [z,w,iw]=ode('stiff',y0,0,0.4,'fex','jex');
  33. z=ode('stiff',z,0.4,4,'fex','jex',w,iw);
  34. if maxi(z-y3(:,2)) > %eps then pause,end
  35.  
  36. [y1,w,iw]=ode(y0,t0,t1(1),'fex');
  37. y2=ode(y0,t1(1),t1(2:nt),'fex',w,iw);
  38. if maxi([y1 y2]-yref) > Leps then pause,end
  39.  
  40. [y1,w,iw]=ode(y0,t0,t1(1),'fex','jex');
  41. y2=ode(y0,t1(1),t1(2:nt),'fex','jex',w,iw);
  42. if maxi([y1 y2]-yref) > Leps then pause,end
  43.  
  44. //   variation of tolerances
  45. atol=[0.001,0.0001,0.001];rtol=atol;
  46. //    externals 
  47. comp(f),comp(j)
  48. //  4. type given , scilab lhs ,jacobian not passed
  49. y4=ode('stiff',y0,t0,t1(1),atol,rtol,f);
  50. if maxi(y4(:,1)-yref(:,1)) > 0.01 then pause,end
  51. //  5. type non given, rhs and scilab jacobian
  52. y5=ode(y0,t0,t1,f,j);
  53. if maxi(y5-yref) > Leps then pause,end
  54. //  6. type given (stiff),rhs and jacobian  by scilab
  55. y6=ode('stiff',y0,t0,t1,0.00001,0.00001,f,j);
  56. if (y6-yref) > 2*0.00001 then pause,end
  57. //  7. matrix rhs, type given(adams),jacobian not passed
  58. // 
  59. a=rand(3,3);ea=exp(a);
  60. deff('[ydot]=f(t,y)','ydot=a*y')
  61. comp(f);
  62. t1=1;y=ode('adams',eye(a),t0,t1,f);
  63. if maxi(ea-y) > Leps then pause,end
  64. //
  65. //   DAE's
  66. //     dy1/dt = -.04*y1 + 1.e4*y2*y3
  67. //     dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
  68. //       0.   = y1 + y2 + y3 - 1.
  69. //  y1(0) = 1.0, y2(0) = y3(0) = 0.
  70. // scilab macros
  71. deff('[r]=resid(t,y,s)','...
  72. r(1)=-0.04*y(1)+1.d4*y(2)*y(3)-s(1),...
  73. r(2)=0.04*y(1)-1.d4*y(2)*y(3)-3.d7*y(2)*y(2)-s(2),...
  74. r(3)=y(1)+y(2)+y(3)-1')
  75. deff('[p]=aplusp(t,y,p)','...
  76. p(1,1)=p(1,1)+1,...
  77. p(2,2)=p(2,2)+1')
  78. deff('[p]=dgbydy(t,y,s)','[-0.04,1.d4*y(3),1.d4*y(2);...
  79. 0.04,-1.d4*y(3)-6.d7*y(2),-1.d4*y(2);...
  80. 1,1,1]')
  81. comp(resid);comp(aplusp);comp(dgbydy);
  82. //        %ODEOPTIONS tests
  83. //       
  84. rand('seed',0);rand('normal');
  85. nx=20;A=rand(nx,nx);A=A-4.5*eye;
  86. deff('y=f(t,x)','y=A*x')
  87. deff('J=j(t,x)','J=A')
  88. x0=ones(nx,1);t0=0;t=[1,2,3,4,5];
  89. nt=size(t,'*');
  90. eAt=exp(A*t(nt));
  91.  
  92. //        Test itask=%ODEOPTIONS(1)
  93.  
  94. //itask=1  --->  usual call (t=vector of instants, solution at all t)
  95. //========================
  96. itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
  97. jacflag=2;ml=-1;mu=-1;
  98. %ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
  99. xf=ode(x0,t0,t,f);   //lsoda
  100. if norm(xf(:,nt)-eAt*x0)>Leps then pause,end
  101. xfj=ode(x0,t0,t,f,j);   //lsoda with jacobian
  102. if norm(xfj(:,nt)-eAt*x0)>Leps then pause,end
  103.  
  104.  
  105. //itask=2;   --->  solution at mesh points  ---> t=tmax
  106. //========================
  107. %ODEOPTIONS(1)=2;tmax=t(5);
  108. xft=ode(x0,t0,tmax,f);
  109. [p,q]=size(xft);
  110. xlast=xft(2:nx+1,q);
  111. if xft(1,q)<tmax then pause,end
  112. if norm(xlast-exp(A*xft(1,q))*x0)>Leps then pause,end
  113.  
  114. //itask=3;   ---> solution at first mesh point beyond t=tmax
  115. %ODEOPTIONS(1)=3;
  116. x3=ode(x0,t0,tmax,f);
  117. if norm(x3(2:nx+1)-xlast,1)>Leps then pause,end 
  118.  
  119. //itask=4; test with %ODEOPTIONS(2)=tcrit
  120. %ODEOPTIONS(1)=4; //---> computation at all t and t>tcrit are not called
  121. tcrit=2.5;%ODEOPTIONS(2)=tcrit;
  122. chk=0;
  123. deff('y=fcrit(t,x)',['if t<=tcrit then'
  124.                       ' y=A*x;'
  125.                       'else'
  126.                       ' y=A*x;chk=resume(1);end'])
  127. x42=ode(x0,t0,t,fcrit);
  128. if chk==1 then pause,end
  129. [p,q]=size(x42);
  130. if norm(x42(:,q)-ode(x0,t0,tcrit,f),1)>Leps then pause,end    
  131.  
  132. //itask=5; test with %ODEOPTIONS(2)=tcrit
  133. %ODEOPTIONS(1)=5;  //---> computation at mesh points and t>tcrit are not called
  134. %ODEOPTIONS(6)=2;  // Estimated jacobian
  135. chk=0;
  136. x52=ode(x0,t0,2.3,fcrit);
  137. if chk==1 then pause,end
  138. [p,q]=size(x52);
  139. if x52(1,q)>tcrit then pause,end
  140.  
  141. //test of %ODEOPTIONS(3:5)=[h0,hmax,hmin]
  142. %ODEOPTIONS(1)=1;
  143. h0=0.0;hmax=0.1;hmin=0.0001;
  144. %ODEOPTIONS(3:5)=[h0,hmax,hmin];
  145. x35=ode(x0,t0,t,f);
  146. if norm(x35-xf,1)>Leps then pause,end
  147.  
  148. //test of %ODEOPTIONS(6)=jacflag
  149. %ODEOPTIONS(6)=1;//Jacobian given
  150. %ODEOPTIONS(3:5)=[0 0 0];
  151. x61=ode('st',x0,t0,t,f,j);   //with Jacobian
  152. if norm(x61-xf,1)>10*Leps then pause,end
  153.  
  154.  
  155. %ODEOPTIONS(6)=0; // jacobian nor called nor estimated
  156. x60=ode('st',x0,t0,t,f,j);   //Jacobian not used (warning)
  157. x60=ode('st',x0,t0,t,f);    //Jacobian not used
  158. if norm(x60-x61,1)>10*Leps then pause,end
  159.  
  160.  
  161. %ODEOPTIONS(6)=1;//Jacobian estimated
  162. x60=ode('st',x0,t0,t,f)  ; 
  163. if norm(x60-x61,1)>10*Leps then pause,end
  164.  
  165. //test of %ODEOPTIONS(6)=jacflag   (adams)
  166. %ODEOPTIONS(6)=1;//with given Jacobian
  167. x60=ode('ad',x0,t0,t,f,j) ;  
  168. if norm(x60-x61,1)>10*Leps then pause,end
  169.  
  170.  
  171. %ODEOPTIONS(6)=0;// jacobian nor called nor estimated
  172. x60=ode('ad',x0,t0,t,f,j);   //Jacobian not used (warning)
  173. x60=ode('ad',x0,t0,t,f);    //Jacobian not used
  174. if norm(x60-x61,1)>10*Leps then pause,end
  175.  
  176. // test lsoda
  177. %ODEOPTIONS(6)=2;// jacobian  estimated
  178. x60=ode(x0,t0,t,f);
  179. if norm(x60-x61,1)>10*Leps then pause,end
  180.  
  181. %ODEOPTIONS(6)=1;//Jacobian estimated
  182. x60=ode(x0,t0,t,f);
  183. if norm(x60-x61,1)>10*Leps then pause,end
  184.  
  185.  
  186. //   Banded Jacobian
  187. itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
  188. //provisional values as default
  189. jacflag=2;ml=-1;mu=-1;
  190.  
  191. %ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
  192. jacflag=2;%ODEOPTIONS(6)=jacflag;   //Banded Jacobian, given
  193. nx=20;A=diag(-[1:nx]);x0=ones(nx,1);
  194. t0=0;t=[1,2,3,4,5];
  195. for k=1:nx-1, A(k,k+1)=1;end
  196. for k=1:nx-2, A(k,k+2)=-2;end
  197. for k=1:nx-1, A(k+1,k)=-1;end
  198. for k=1:nx-2, A(k+2,k)=2;end
  199. for k=1:nx-3, A(k+3,k)=-3;end
  200. deff('xd=f(t,x)','xd=A*x')
  201. ml=3;mu=2;
  202. %ODEOPTIONS(11:12)=[ml,mu];
  203. for i=1:nx;
  204.     for j=1:nx;
  205. if A(i,j)<>0 then J(i-j+mu+1,j)=A(i,j);end
  206. end;end;
  207. // J is a ml+mu+1 x ny matrix.
  208. // Column 1 of J is made of mu zeros followed by df1/dx1, df2/dx1, df3/dx1,
  209. // i.e. 1 + ml possibly nonzeros entries.
  210. // Column 2 of J is made of mu-1 zeros followed by df1/dx2, df2/dx2, ...
  211. // etc...
  212. %ODEOPTIONS(6)=1;%ODEOPTIONS(11:12)=[-1,-1];
  213. deff('jj=j1(t,x)','jj=A')
  214. xnotband=ode('st',x0,t0,t,f,j1);
  215.  
  216.  
  217. %ODEOPTIONS(6)=4;//banded jacobian external given
  218. %ODEOPTIONS(11:12)=[3,2];
  219. deff('jj=j(t,x)','jj=J')
  220. xband=ode('st',x0,t0,t,f,j);
  221. if norm(xnotband-xband,1)>Leps then pause,end
  222.  
  223. %ODEOPTIONS(6)=5;//banded jacobian evaluated
  224. %ODEOPTIONS(11:12)=[3,2];
  225. deff('jj=j(t,x)','jj=J')
  226. xband=ode('st',x0,t0,t,f,j);
  227. if norm(xnotband-xband,1)>Leps then pause,end
  228.  
  229.  
  230. //            Test of %ODEOPTIONS(7)  
  231. //%ODEOPTIONS(7)=mxstep  ---> maximum number od steps allowed
  232. itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
  233. //provisional values as default
  234. jacflag=2;ml=-1;mu=-1;
  235. %ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
  236. %ODEOPTIONS(7)=10;
  237. //ode(x0,t0,t,f);  // ---> Non convergence
  238.  
  239. //            Test of %ODEOPTIONS(8:9)  
  240. //%ODEOPTIONS(8:9)=[maxordn,maxords] ---> maximum order for nonstiff and stiff
  241. itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
  242. //provisional values as default
  243. jacflag=2;ml=-1;mu=-1;
  244. %ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
  245. %ODEOPTIONS(8:9)=[0,0]; //---> default values
  246. wref=ode(x0,t0,t,f); //just for computing reference
  247.  
  248. %ODEOPTIONS(8:9)=[4,3];
  249. ww=ode(x0,t0,t,f);norm(wref-ww,1);
  250.  
  251. %ODEOPTIONS(8:9)=[12,5];
  252. if norm(wref-ode(x0,t0,t,f),1)>Leps then pause,end
  253.  
  254. //using stiff method
  255.  
  256. %ODEOPTIONS(9)=0;
  257. wref=ode('st',x0,t0,t,f);
  258. %ODEOPTIONS(9)=5;
  259. if norm(wref-ode('st',x0,t0,t,f),1) >Leps then pause,end //=0
  260. %ODEOPTIONS(9)=4;
  261. if norm(wref-ode('st',x0,t0,t,f),1)  >Leps then pause,end  //small
  262.  
  263.  
  264. //using nonstiff method
  265.  
  266. %ODEOPTIONS(8)=0;
  267. wref=ode('ad',x0,t0,t,f);
  268. %ODEOPTIONS(8)=12;
  269. if norm(wref-ode('ad',x0,t0,t,f),1) >Leps then pause,end  //=0
  270. %ODEOPTIONS(8)=5;
  271. if norm(wref-ode('ad',x0,t0,t,f),1) >Leps then pause,end   //small
  272.  
  273. //mixed 
  274. %ODEOPTIONS(8:9)=[5,12];
  275. wref=ode(x0,t0,t,f);
  276. %ODEOPTIONS(8:9)=[4,10];
  277. if norm(ode(x0,t0,t,f)-wref,1)>Leps then pause,end   //small
  278.  
  279.  
  280. A=diag([-10,-0.01,-1]);
  281. deff('uu=u(t)','uu=sin(t)');
  282. B=rand(3,1);
  283. deff('y=f(t,x)','y=A*x+B*u(t)')
  284. %ODEOPTIONS(1)=2;
  285. yy1=ode('stiff',[1;1;1],0,1,f);
  286. yy2=ode('stiff',[1;1;1],0,2,f);
  287. %ODEOPTIONS(1)=3;
  288. yy1=ode('stiff',[1;1;1],0,1,f);
  289. yy2=ode('stiff',[1;1;1],0,2,f);
  290.  
  291. clear %ODEOPTIONS;
  292. rand('seed',0);rand('normal');
  293. nx=20;A=rand(nx,nx);A=A-4.5*eye;
  294. deff('y=f(t,x)','y=A*x')
  295. deff('J=j(t,x)','J=A')
  296. //%ODEOPTIONS(1)=1;
  297. y2=ode('stiff',ones(nx,1),0,2,f,j);
  298. [y1,w,iw]=ode('stiff',ones(nx,1),0,1,f,j);
  299. y2p=ode('stiff',y1,1,2,f,j,w,iw);
  300. y12=ode('stiff',ones(nx,1),0,[1,2],f,j);
  301. norm(y12(:,2)-y2p);
  302. yaf=ode('adams',ones(nx,1),0,2,f,j);
  303. yaj=ode('adams',ones(nx,1),0,2,f,j);
  304. ysf=ode('stiff',ones(nx,1),0,2,f,j);
  305. ysj=ode('stiff',ones(nx,1),0,2,f,j);
  306.  
  307.  
  308.